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Abstract. Various types of expansions in series of Cheby- 
shev-Hcrmite polynomials currently used in astrophysics 
for weakly non-normal distributions are compared, namely 
the Gram-Charlier, Gauss-Hermite and Edgeworth expan- 
sions. It is shown that the Gram-Charlier series is most 
suspect because of its poor convergence properties. The 
Gauss-Hermite expansion is better but it has no intrinsic 
measure of accuracy. The best results are achieved with 
the asymptotic Edgeworth expansion. We draw attention 
to the form of this expansion found by Petrov for arbitrary 
order of the asymptotic parameter and present a simple al- 
gorithm realizing Petrov's prescription for the Edgeworth 
expansion. The results are illustrated by examples similar 
to the problems arising when fitting spectral line profiles 
of galaxies, supernovae, or other stars, and for the case of 
approximating the probability distribution of peculiar ve- 
locities in the cosmic string model of structure formation. 
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1. Introduction 

The normal, or Gaussian, distribution plays a prominent 
role in statistical problems in various fields of astrophysics 
and general physics. This is quite natural, since the sums 
of random variables tend to a normal distribution when 
the quite general conditions of the central limit theorem 
are satisfied. In many applications, to extract useful in- 
formation on the underlying physical processes, it is more 
interesting to measure the deviations of a probability den- 
sity function (hereafter PDF) from the normal distribu- 
tion than to prove that it is close to the Gaussian one. This 
has been done for example in the work on peculiar veloci- 
ties and cosmic microwave background anisotropies in var- 
ious cosmologi cal m odels (Scherrer Sz Ber tschinger 1991 
Kofman et al 



1994; Moessner et al. 1994; Juszkiewicz et 



al. 1995; Bernardeau & Kofman 1995. Amendola 1994 



Colombi 1994; Moessner 1995; Ferreira et al. 1997; Gaz- 



tanaga et al. 1997), in analyzing the velocity distributions 
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and fine structure in ellip tical galaxies ( Rix fc White 
van d er Marel & Franx |l993| ; Gerhard |1993| ; Heyl et al. 
1994 ) , and in st udies of large Reynolds number turbulence 



(Tabeling et al. 1996 ). In this paper we present a unified 
approach to the formalism used in those applications, il- 
lustrating it by examples of similar problems arising in 
cosmology and in the theory of supernova line spectra. 

The first investigation of slightly non-Gaussian distri- 
butions was undertaken by Chebyshev in the middle of the 
19th century, who studied in detail a family of orthogonal 
polynomials which form a natural basis for the expansions 
of these distributions. A few years later the same poly- 
nomials were also investigated by Hcrmitc and they are 
called Chebyshev-Hermite or simply Hermite polynomials 
now (their definition was first given by Laplace, see e.g. 
Encyclopaedia of Mathematics 1988). 



There are several forms of expansions using Hermite 
polynomials, namely the Gram-Charlier, Gauss-Hermite 
and Edgeworth expansions. We introduce the notation in 
Sect. H and use the simple example of the \ 2 distribution 
for various degrees of freedom to illustrate the properties 
of Gram-Charlier expansions in Sect. ||. In subsequent sec- 
tions the x 2 distribution is used for testing the other two 
expansions. We note in Sect. || that the Gram-Charlier 
series is just a Fourier expansion which diverges in many 
situations of practical interest, whereas the Gauss-Hermite 
series has much better convergence properties. However, 
we point out that another series, the so called Edgeworth 
expansion, is more useful in many applications even if it is 
divergent, since, first, it is directly connected to the mo- 
ments and cumulants of a PDF (the property which is 
lost in the Gauss-Hermite series) and, second, it is a true 
asymptotic expansion, so that the error of the approxima- 
tion is controlled. 



ner 



Some applications (e.g. Bernardeau 1992,1994; Moess- 
1995] ) involve cumulants of higher order. These require 



the use of a correspondingly higher order Edgeworth se- 
ries, but only the first few terms of the series are given 



in standard references (Cramer 1957; Abramowitz & Ste 
Juszkiewicz et 



gun 
man 



Bernardeau 



Kof- 



1972 ; Juszkiewicz et al. 1995 

1995 ). In Sect. || we popularize a deri v ation of the 
Edgeworth expansion due to Petrov (|962L |l972| , |l987| ) 
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for an arbitrary order of the asymptotic parameter, and 
we present a slightly simpler and more straightforward 
way of obtaining his result. The formula found by Petrov 
(see Sect. ^ below) requires a summation over indices with 
non-trivial combinatorics which hindered its direct appli- 
cation. We find a simple algorithm realizing Petrov's pre- 
scription for any order of the Edgeworth expansion and 
this algorithm is easily coded, e.g. with standard Fortran, 
eliminating the need for symbolic packages. 

We find that the Edgeworth-Petrov expansion is in- 
deed very efficient and reliable. In Sect. ^| we apply the 
formalism to the problem of peculiar velocities resulting 
from cosmic strings studied by Moessner ( |1995 ) and we 
show how this technique allows one to reliably extract de- 
viations from Gaussianity, even when they are tiny. 



2. Background and notation 

Let F(x) be the cumulative probability distribution of a 
random variable X. Then the mean value for the random 
variable g(X) is the expectation value 



Eg(X) ee (g(X)) = / g(x)dF(x) 



(1) 



The PDF is p(x) = dF(x)/dx. The distribution F(x) 
is not necessarily smooth, so it can happen that p(x) 
is nonexistent at certain points. Nevertheless, the mean 
Eg(X) is defined as long as the integral in (Q) exists. Fol- 
lowing the definition (m) , the fc-th order moment of X is 



Ctk 



EX K = 



x k dF{x) . 



(2) 



Thus, the mean of X is m = ot\ = EX and its dispersion 
is 



/oc 
(x - m) 2 dF{x) . 
-oc 

We denote the cumulative normal distribution by 

f x 1 f x 

P(x) ee / Z(t)dt = —= I cxp(-t 2 /2)dt , 

J — oo \2lT J — qo 

so its PDF is the Gaussian function 
exp(— x 2 /2) 



Z(x) 



/27T 



(3) 



(4) 



(5) 



We also need to recall some definitions for sets of or- 
thogonal polynomials. Any two polynomials P n (x) and 
P m (x) of degrees n ^ m are orthogonal on the real axis 
with respect to the weight function w(x) if 



w(x)P n (x)P m (x)=0 



(6) 



We follow the notation of Abramowitz & Stegun ( 1972 ) 
where possible, so He n {x) denotes the polynomial with 



weight function w(x) = exp(— x 2 /2) cx Z(x). According 
to Rodrigues' formula, 



He n (x) = (-ire* 2 /**-* 2 /! 



dx n 



(7) 



We will refer to He n (x) as Chebyshev-Hermite polynomi- 
als following Kendall ( |1952| ). The ones with the weight 

w(x) = exp(— x 2 ) oc Z 2 (x) 1 



H n (x) = (-l)«e» ^e"* , 



(8) 



will be called Hermite polynomials. 

The expansions in Chebyshev-Hermite and Hermite 
polynomials are used in many applications in astrophysics. 
For a PDF p(x) which is nearly Gaussian, it seems natural 
to use the expansion 



p{x) 



^ d n Z(x) 
^ Cn dx n 

n=0 



(9) 



where the coefficients c„ measure the deviations of p(x) 
from Z(x). From the definitions (||) and (^) it follows that 



d n Z(x) 

dx n 
so that 

p(x) - ' 



= {-l) n He n (x)Z(x) , 



(-l) n c n He n {x)Z{x) , 



with 



(-1)" 



p{t)He n (t)dt 



(10) 



(11) 



(12) 



This is the well-known Gram-Charlier series (of type A, 
see e.g. Cramer 1957; Kendall 1952 and references therein 
to the original work). Since He n {x) is a polynomial, we 
see that the coefficient c„ in ([l2]) is a linear combination 
of the moments of the random variable X with PDF 
p{x). The combination is easily found by using the explicit 
expression for He n (x), which we derive in[H, namely 



He n {x) 



= n.\ 



(-1) 



^-^ k\(n 



2k)\ 2 fc 



(13) 



The Gram-Charlier series was used in cosmological ap- 
plications in the paper by Scherrer & Bertschinger (1991), 
but note that it is incorrectly called Edgeworth expansion 
there. The Gram-Charlier series is merely a kind of Fourier 
expansion in the set of polynomials He n (x). This expan- 
sion has poor convergence properties (Cramer 1957). Very 
often, for realistic cases, it diverges violently. We consider 
such an example in the next section. The Edgeworth se- 
ries, discussed in detail in Sect. ||, is a true asymptotic 
expansion of the PDF, which allows one to control its ac- 
curacy. 



S. Blinnikov & R. Moessner: Expansions for nearly Gaussian distributions 



3 



3. An example based on the x 2 distribution 

A good example for illustrating the fast divergence of the 
Gram-Charlier series is given by its application to the xt 
distribution with v degrees of freedom, since the moments 
of this distribution are known analytically and its PDF 
tends to the Gaussian one for large v. If Xi, X2, ■ X„, are 
independent, normally distributed random variables with 
zero expectation and unit dispersion, then the variable 
X 2 = xl = E:=i X? has the PDF 



0.4 



p(x 2 



(x 2 r /2-l exp( _ x 2 /2) 

2^/ 2 r(i//2) 



x 2 >o 



(14) 



0.2 



The expectation value of x 2 is v an d its dispersion is 2v. 
If x = (x 2 — v) /v2^j then x has zero expectation and unit 
dispersion and its PDF p(x) asymptotically tends to the 
Gaussian distribution Z(x). Transforming p(x 2 ) from ( |l4|) 
to p(x), we obtain for x > —\J v/2: 



Yj & N = 2 in Gram-Charlier 




p(x) 



<2v 



(V2vx + v)"' 2 - 1 exp ~(s/2vx + v)/2 

2^/ 2 r(t//2) 



(15) 



0.4 



0.2 



The x 2 distribution was employed by Matsubara & Yoko- 
yama ( 1996[ ) for a representation of the cosmological den- 
sity field (see also Luo 1995| for an application of the x 2 dis- 
tribution in the context of cosmic microwave background 
temperature anisotropies). It can also serve as a crude 
model for approximating the profiles of spectral lines in 
moving media. If the lines are nearly Gaussian in the mat- 
ter at rest, then the motion with a high velocity gradi- 
ent, like that in supernova envelopes, leads to distortions 
which can be approximated by the xl PDF (with v = 2 
for highly non-Gaussian profiles, see e.g. Blinnikov 1996, 
and references therein). 

The comparison of the Gram-Charlier approximations 
of p(x) with the exact results is presented in Figs. |l] and 
|^ for an increasing number of terms in the expansion. It 
is clear that the series quickly becomes inaccurate with a 
larger number of terms included. 

4. Fourier expansions 

In order better to understand the poor convergence prop- Fig. 1. The normalized \ 2 PDF ( |l5| ) for v = 5 (dashed 

erties of the Gram-Charlier series, let us first discuss how line), and its Gram-Charlier approximations with 2 and 6 

it is rel ated t o the Fo urier expansion. A Fourier expan sion terms in the expansion (solid line) 
(Szego |1978[ Suetin |l97S| ; Nikiforov & Uvarov 




for 

any function f(x) in the set of orthogonal polynomials P n 
is given by 



/(*) 



n=0 



Q>n-Pn(x*) , 



with 



— / w(t)f(t)P n (t)dt 

^"n J —00 



(16) 



(17) 



Here h n is the squared norm 

/oo 
w{t)Pl{t)dt , 
-00 

and 

h n = V2nn\ for He n (x) , 
h n = y/^2 n nl for H n (x) . 



(18) 

(19) 
(20) 
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"ft 



10 



-20 - 



= 12 in Gram-Charlier 




2xl0 13 



-2xl0 13 - 




Fig. 2. The same as in Fig. [j] but for 12 and 36 terms in 
the Gram-Charlier expansion. 



Now we can see that the Gram-Charlier series (|llj) is just 
the Fourier expansion © of f(x) = p(x)/Z(x) in the set 
of Chebyshcv-Hcrmitc polynomials with c„ = (— l) n a n . 

The properties of the Gram-Charlier approximations 
of p(x) in Figs. [I] and ^| are to be considered in the gen- 
eral context of the convergence of Fourier expansions. The 
source of the divergence lies in the sensitivity of the Gram- 
Charlier series to the behavior of p(x) at infinity - the lat- 
ter must fall to zero f aster than exp( — x 2 / A) for the series 



too restrictive for practical applications. Our example of 
the x 2 distribution in ([l^), with its exponential behavior 
at infinity, clearly demonstrates this. 

The Fourier expansion of p(x)jZ{x) in another set 
of Hermite polynomials H n (x) (||) (not in Chebyshev- 
Hermite polynomials He n (^), as for the Gram-Charlier 
series) is sometimes used: 



P( x ) ~ X! a nH n {x)Z(x) , 



with 



71 = 



Z{t)p(t)H n (t)dt 



(21) 



(22) 



This series is often called the Gauss-Hermite expansion 
(see e.g. its application to spectral lines of galaxies in van 
der Marel & Franx 1993). Examples in Figs. || and || show 
its better convergence. 



Van der Marel & Franx (1993) use a theorem due to 
Myller-Lebedeff on the convergence of the Gauss-Hcrmitc 
expansion: it converges when x 3 p(x) — > for any p{x) 
with finite and continuous second derivative. Actually, the 
conditions sufficient for the convergence are better: if p(t) 
obeys the Lipschitz condition 



\p(t)-p(x)\ <M\t- 



M = const, < a < 1 , (23) 



in a vicinity of x and 

b(i)|(l + |i| 3/2 )*<(X3, 



(24) 



then the G auss-H ermite series converges to p{x) at x (see 
e.g. Suetin 1979 ). These weaker conditions imply that the 
class of PDFs with convergent Gauss-Hermite expansions 
is much wider than suggested by the Myll er-Lebedeff theo- 
rem cited in van der Marel & Franx 1993 ). But the simple 
relation between the coefficients in the expansion and the 
moments of the PDF typical for the Gram-Charlier se- 
ries is now lost [cf. Eq. 



to converge (Cramer 1957; Kendall 1952). This is often 



L2J), (|13|) and (|22|)]. It might be 
not important in many practical applications when the 
moments cannot be accurately determined from observa- 
tions, but it is very important for a theoretical work based 
on the analysis of the moments. 

Our Figs. H and [| show that the \ 2 PDF is well-suited 
for approximation by a Gauss-Hermite series. However, 
one should be cautious about the accuracy of the compu- 
tation of the Fourier coefficients for higher order terms. 
We have not encountered the problem of numerical errors 
in our Gram-Charlier example, since there all coefficients 
can be calculated analytically. Yet in general care must 
be taken of the computational accuracy to avoid spurious 
numerical divergence of a series which converges theoret- 
ically. 
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Fig. 3. The normalized x 2 PDF ( |l5| ) for v = 5 (dashed Fig. 4. The same as in Fig. || but for 12 and 36 terms in 
line) , and its Gauss- Hcrmitc approximations with 2 and 6 the Gauss- Hermite expansion 
terms in the expansion (solid line). 



5. The Edgeworth asymptotic expansion 

A random variable X can be normalized to unit variance 
by dividing by its standard deviation a. The Edgeworth 
expansion is a true asymptotic expansion of the PDF of 
the normalized variable X/a in powers of the parameter 
a, whereas the Gram-Charlier series is not. This difference 
between the Gram-Charlier series and the Edgeworth ex- 
pansion was pointed out by Juszkiewicz et al.(1995), and 



in independent work by Bernardeau & Kofman (1995), 
followed by Amendola fll994p and Colombi fll994| ). 

Juszkiewicz et al. ( 1995| ) and Bernardeau & Kofman 



( 1995 ) use 2 or 3 ter ms of the Edgeworth expansion de- 
rived e.g. in Cramer ( 1957 ). We note that the full explicit 
expansion for arbitrary order s was obtained already in 
1962 by Petrov ( |l962| ), see also Petrov flL972j |1987[ ). 

Petrov derived a powerful generalization of the Edge- 
worth expansion for a sum of random variables. In this 
section we give a simplified derivation of the Edgeworth se- 
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Fig. 5. The normalized x 2 PDF ( |Hf ) for v = 2 (dashed Fig. 6. The normalized \ 2 PDF (|l5|) for v = 5 (dashed 
line), and its Edgeworth-Petrov approximations with 4 line), and its Edgeworth-Petrov approximations with 2 
and 12 terms in the expansion (solid line). and 4 terms in the expansion (solid) 



ries, following Petrov ( 1972 ). This derivation, for an arbi- 
trary order, is somewhat simpler than the derivation given 
for example by Bernardeau & Kofman (1995) for the third 
order only of the Edgeworth expansion. 

The characteristic function Q(t) of a random variable 
X is the expectation E exp(itX) as a function of t, 



*(t) = 



c dF(x) 



(25) 



that is the Fourier transform of p(x) if the probability den- 
sity p(x) = dF(x)/dx exists. The definition ( |25| ) implies 
that if the moment (0) of X exists, 



$W(0) = i k a k . 

Hence the Taylor series for <!>(i) is given by 



*(*)~i+E^( tt ) 



k=l 



(26) 



(27) 
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A similar series for ln<I>(i), Let us write the exponential function in (Q) as a formal 

series in powers of a, 



involves cumulants (semi-invariants) « n defined by 

where the coefficient of the power s is a function V s (it). 



1 



7" 



f-ln*(t) . (29) Now, using 5 (a;) = E" i{5 r+ 2^r +2 a;7^ + 2)!} and/ = 



t=o 



In ^ we prove a fundamental lemma of calculus for the 
n— th derivative of a composite function fog(x) = f(g(x)), V s (it) = ———f(g(x))\ x =o 
which reads s ' 



exp in (|30|), we find that 
I d s 



S m+2 (it) m + 2 \ k " 

/(<?(*)) = "^HWl (m+2)! 



d n 

dx" J (/■•„,! J 

n! 



1/1 \ 

/' r H2/)lt/=9(x) ( — ff^M ) j (30) where the summation extends again over all non-negative 

m —i / integers {k m } satisfying (|3l"|). Thus the function P s is just 

a polynomial. 

where r = fa + k 2 + ... + k n and the set {£;,„} consists Suppose that the probability density p(x) of a ran- 

of all non-negative integer solutions of the Diophantine dom variable X exists. Then the PDF for X/a is q(x) = 

equation ap(ax), and it is the inverse Fourier transform of the char- 



fa + 2k 2 + ... + nk n = n . (31) 



acteristic function <p: 



1 r°° _. 

Using @ we derive a useful relation (Petrov [1987D = — / e ltx ip{t)dt . (38) 

between the cumulants k„ and the moments of a PDF v 



m 



^.^(-inr-Dinri© 1 "' < 32 > 



If <I>(t) is the Fourier transform of a function p(x), then 
(— it) n $(t) is the transform of the n-th derivative of p(x), 



m=i— ■ ±_ p{x) = _J e -^(-it) n <l>(t)dt. (39) 
Here summation extends over all non-negative integers 

{k m } satisfying @ and r = fci + fc 2 + ... + fc„. We describe The Fourier transform of the Gaussian distributio n Z(x ) 

a simple algorithm for obtaining all solutions of Eq. @ in (§) is exp(-i 2 /2) , see e.g. Bateman & Erdelyi ( |1954{ ). 

i n R, Therefore each (it)™, multiplied by exp(— i 2 /2) in the ex- 

Now we are ready to begin with the derivation of the pansjon of <p (see Eqs. |§ to |37|), generates (according to 

Edgeworth expansion. Consider a random variable X with Eq. |38|) the n-th derivative of 
EX = (this can always be achieved by an appropri- 

ate choice of origin), and let X have dispersion a 2 . If (_i)™_l_^(V) = / e - ttx (it) n exp(~t 2 /2)dt , (40) 
X has the characteristic function $(i), then the normal- ./-oo 
ized random variable X/a has the characteristic function 

p(t) = $(t/a). Therefore we have from Eqs. (§) and © m the corresponding expansion for q{x), 
that oo 



oo 

ln^)=ln$(t/a)~^^-(^«. (33) 




Here the sum starts at n = 2 because EX = 0. Moreover, | V (m + 2)! dx m+2 

since «2 = cr 2 (see Eq. (|32"|)) we obtain 

{Here the set |fc m j in the sum consists of all non-negative 

g bnO_ ( . f)n ^ (34) integer solutions of the equation 

™ =3 fe x + 2fe 3 + ... + sfc a =s. (42) 

with 



Using ( |l0| ) and r = fci + fc 2 + •■■ + k s we can rewrite (jiH ) 
S n = n n /a 2n ~ 2 . (35) in terms of the Chebyshev-Hermite polynomials: 
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X% & N=12 in Edgeworth 



0.5 



x 

"ft 



-0.5 




0.4 



X 



0.2 



X20 & N = 2 in Edgeworth 




_j 1 1 1 1 



-2 



Fig. 7. The normalized x 2 PDF (|Ht ) for v = 5 (dashed 
line) and its Edgeworth-Petrov approximations with 12 
terms in the expansion (solid) 



q(x) = ap{ax) = Z(x) I 1 + ^ a s 



(43) 



This is the Edgewort h expansion for arbitrary order 
s. See Petrov ( 1972 , 1987 ) for a more general form of the 
expansion (for non-smooth cumulative distribution func- 
tions F(x) and for a sum of random variables) and for the 
proof that the series ( f43|) is a symptotic (see also the clas- 
sical references Cramer 1957 and Feller [1966 ). This means 
that if the first ./V terms are retained in the sum over s, 
then the difference between q(x) and the partial sum is 
of a lower order than the N-th term in the sum (Erdelyi 



1956; Evgrafov 1961). Convergence plays no role in the 



definition of the asymptotic series. 

Strictly speaking, Petrov (1972) proves the asymptotic 
theorems for sums of v independent random variables only 
when a ~ and not for any a , which we used 

in our derivation. But in all practical applications where 
nearly Gaussian PDFs occur (and in all applications that 
we consider in the present work), those PDFs basically 
are the sums of random variables, and the proofs of the 
asymptotic theorems are relevant. In the next section we 
show how the theory works in practice. 

Figures || - 1| show some examples of the Edgeworth ex- 
pansion for the x 2 distribution. It is clear that for strongly 



0.4 



0.2 




Fig. 8. The normalized X 2 PDF @ for v = 20 (dashed 
line), and its Edgeworth-Petrov approximations with 2 
and 4 terms in the expansion (solid) 



non-Gaussian cases, like xt f° r v = 2, it has a very small 
domain of applicability in practical cases since it diverges 
like the Gram-Charlier series for a large number of terms 
(Fig. [sj> . But already in this case one can check that the 
order of the last term retained gives the order of the error 
correctly, and one can truncate the expansion when the 
last term becomes unacceptably large. For nearly Gaus- 
sian distributions the situation is much better: compare 
the cases for v = 5 and v = 20 in Figs. ^|, and |[ 
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Fig. 9. Edgeworth expansion up to 10th order for the PDF Fig. 10. Relative deviation E^r from a normal distribution 
of peculiar velocities from cosmic strings, within an ana- of the Edgeworth expansion up to A~th order of the PDF 
lytical model for the string network, for two values of the of peculiar velocities from cosmic strings, for n 8 = 1. Also 

shown is the error tjy of this deviation associated with the 
expansion. 



number of strings per Hubble volume (Moessner 1995) 



6. Peculiar velocities from cosmic strings 



For simplicity, let us write Eq. schematically as 



As an example we consider the probability distribution 
of peculiar velocities within the cosmic string model of 
structure formation. Cosmic strings are one-dimensional 
topological defects possibly formed in a phase transition 
in the early Universe (Brandenberger 1994; Hindmarsh & 
Kibble |1995[ Vilenkin & Shellard [L994Q. After the time 



q(x) = Z(x){l 



of formation, the string network quickly evolves to a scal- 
ing solution with a constant number n s of strings passing 
through a Hubble volume in one expansion time. 

Since individual topological defects give rise to velocity 
perturbations of the dark matter through which they move 
up to a distance of a Hubble radius, one expects a non- 
Gaussian result, in contrast to inflationary theories which 
predict a Gaussian PDF. Therefore departures from a nor- 
mal distribution may be a way to distinguish between the 
two main classes of theories of structure formation, infla- 
tion and topological defects. However, since many strings 
present between the time of equal matter and radiation 
contribute to the perturbations, a nearly Gaussian PDF 
can result due to the central limit theorem. 

In order to estimate the deviation of the PDF from the 
normal distribution we use Petrov's formula (|4^) for the 
Edgeworth series in this section. We calculate the cumu- 
lants up to 12th order in the analytical mod el for the cos- 
mic string network presented in Moessner (1995), where 
the cumulants are given up to 8th order. 



s=l 



(44) 



i.e. denote the sth term in the sum by t s . Let us further 
denote the sum up to s — N by Ejv = Yl s =i ^ n Fig- H 
we show the Edgeworth expansion of the PDF of peculiar 
velocities up to 10th order, for the case of n s =l string per 
Hubble volume. Expanding up to A^th order, the relative 
deviation of the PDF from a normal distribution is given 

by, 



1 = E 



(45) 



Z{x) 

It is only significant if the error of the asymptotic expan- 
sion, which is of the order of the last term included, is 
smaller than this deviation. In Fig. [lO] we show the rel- 
ative deviation E^v and the error ijy associated with it. 
The wiggles or cusps in the graphs appear at zeros of Ejv 
and 1 7v which are both oscillating, changing their sign, and 
we plot the logarithms of the absolute values. So only the 
maxima of the curves give a true indication of the devi- 
ations and errors. For N — 10 the error is clearly below 
the deviation and thus significant. For N = 4, the error 
is still below the deviation for most values of x, but it is 
not as clear, especially since the error is not exactly equal 
to tff but of that order only. In Fig. [n] we show the rel- 
ative error t^/{l + Ejy) in the expansion of q{x)/Z(x). 
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Fig. 11. Relative error ^/(l + Ejv) in the Edgeworth 
expansion of q(x)/Z(x) for the PDF of peculiar velocities 
from cosmic strings, for two values of N and n s . 

It is smaller for an expansion to higher order for a given 
number of strings per Hubble volume. The error decreases 
more strongly with N for larger n s . The relative error is 
also smaller, at fixed N, for a larger n Sl in which case the 
PDF is closer to a normal distribution. 

It is interesting t o com pare these results with the gen- 
eral theory. Petrov ( 1972 ) proves that under certain con- 
ditions (which are fulfilled in most physically important 
cases) 



q(x)=Z(x){l+^t s } + o(a") 



(46) 



uniformly for — oo < x < +oo, when a ~ X/v 1 / 2 and the 
PDF q(x) is for a sum of v random variables. One should 
remember that each t s is of the order of a s in (|4^) . In our 
case v is just n s , so the error of the truncated Edgeworth 
series scales as ~ 1/n^ 2 . From Fig. [ll] we can see that 
the error for the case of n s — 10 strings is indeed N/2 
orders less than for the case of n s = 1 string, i.e. 2 orders 
for TV = 4 and 5 orders for N — 10 terms in the expansion. 
This is an illustration of the theory developed by Petrov 



1972,1987) 



7. Conclusions 

We have shown that the Gram-Charlier series has a lim- 
ited domain of applicability for nearly normal distribu- 
tions because of its rather poor convergence properties. 
The Gauss-Hermite expansion can give good results in 



problems like fitting profiles of spectral lines of galaxies, 
supernovae, or ordinary stars. In advanced calculations of 



stellar atmospheres (e.g. Hauschildt et al. 1997; Hubeny 
& Lanz 1995| ) the profiles of thousands or even millions 
of lines must be integrated for up to hundreds of Doppler 
widths, and the Gauss-Hermite expansion can perhaps be 
useful for saving information of the line profiles in an eco- 
nomical way. But since it has no intrinsic measure of accu- 
racy, the number of terms needed in the expansion must 
be examined carefully for each individual problem. 

For situations where the estimate of a deviation of a 
PDF from a Gaussian one is needed, the asymptotic Edge- 
worth expansion is indispensable, and for high order mo- 
ments the form of this expansion found by Petrov is neces- 
sary. We found a workable algorithm for Petrov's formula 
of the Edgeworth expansion and applied it to several ex- 
amples. The source codes used in this work are available 
on request from the authors. 
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Appendices 
Appendix A: Lemma 

In Sect. |, the relation (|30| ) for the n— th derivative of a 
composite function fog{x) = f(g(x)) is used for the deriva- 
tion of the Edgeworth asymptotic expansion. Here a sim- 
plified derivation of Eq.© is given. Petrov ( |1972| , |1987D 
suggests a proof by induction. We note that this deriva- 
tion is more transparent if one simply considers the Taylor 
expansion for fog expressed in terms of the Taylor expan- 
sions of / and g - see Bourbaki (1958). We have 

/( „) 



M = f(yo) 

and 

g(x) = g(x )- 



t 
1! 



g 1 
IT 



Ay- 



r 

2! 



AV 



-AV 



2! 



AV 



9 



(»0 



(AT) 



.(A.2) 



Truncating the expansions at some n and m we find that 

f°g{x) = f(g(x )) 



g 



(m) 



-A r 



+ • 



(A.3) 



/( n) 



1! 



Ax- 



(m) 



2! 



-A' 
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On the other hand, we can write down the Taylor series 
for the composite function, 



f°g(x) = f°g(x ) 

I (/ ° g) '-Ax I {f ° 9) "^ x I ifo9)iS) 
+ 1! + 2! + s! 

Now using the polynomial theorem, 

(xi + x 2 H h x m ) r 

EH fei fe2 



•(A.4) 



{fc m } 



(A.5) 



where summation extends over all sets of non-negative in- 
tegers {k m } satisfying r = k± + k 2 + . . . + k s , and comparing 



the terms A s x with equal s in ( A. 4 ) and (A4), we obtain 
d s 



dx 



■f(g(x)) 



ifi 



{k m } 



h I 



(A.6) 



This is the relation ( pQ) which we sought. Here the set 
{k m } consists of non-negative solutions of the Diophantine 
equation 



ki + 2k 2 + ... + sk s = s , 
and r — ki + k 2 + ... + k s 



(A.7) 



Appendix B: Applications of the Lemma 

If we apply lemma ( ^0| ) to Chebyshev-Hermite polynomi- 
als in (0), we g(x) = —x 2 /2 for / = exp, so that only the 
terms with m = 1 and m = 2 are non-zero in the product 
in (|30|), and we only need non-negative integers {fci, k 2 } as 
the solutions for k\ + 2k 2 = n. Thus for each k = k 2 run- 
ning from to [n/2] (entier of n/2) we have ki = n — 2k 
and r = n — k. Finally, we have from (pQ) that 



d n 
dx™ 



exp{-x 2 /2) 



1! 



(-x) 



[n/2] 
k=Q 

—2k l , l 



1 



k\ \ 2\ 



{n-2k)\ 
k 



r(-i) 



(B.l) 



and the explicit expression ( |l3| ) follows immediately from 
Rodrigues' formula (|^). 

Among other consequences of (|30j) is the relation ( |32"| ) 
between cumulants K n and moments ctk of a PDF. From 
the definition ( p9[ ) we obtain this relation by simply ap- 
plying ( |30| ) to the case of / = In and g = $. Since 

f (r) (y)\y= g ( t) = (-iy-Hr-iy-/$ r \t= - (-ly-v-i)!, 



we find that 

K " = ^dF ln$|t = = 



{k m } 

Thus, from (|2 



(B.2) 



-^(-n-i»nSS)' 



(B.3) 



{k m } 



which is equivalent to (|32|). Here the sum extends over 
all non- negative integers {k m } satisfying ( |3l| ) and r = 
ki + k 2 + ... + k n . 



Appendix C: Algorithm for computing indices k m 

To find all the solutions of Eq. (|l]) it is desirable to order 
all sets of non- negative integers {k m } satisfying it. We first 
rewrite (|lj) as 



2k 2 + k\ = n . 



(C.l) 



It is natural to establish the ordering of the sets Sj = {k m } 
satisfying Eq. ([n]) according to the order of numbers in 
their decimal representation. For n = 3, say, we have 3 
non-negative solutions 



Si 


= {k 3 


= 0,fc 2 


= Q,ki 


= 3} 


s 2 


= {k 3 


= 0,k 2 


= l,fci 


= 1} 


s 3 


= {k 3 


= l,k 2 


= Q,h 


= 0} 



which can simply be written as 

Si = 003, = Oil, #3 = 100 , 

and we say that Si < S 2 < S3 since 3 < 11 < 100. For 
n > 10, when the base 10 is no longer convenient, the sets 
of solutions can be ordered according to the order of the 
integer numbers 



k n r n 1 + ... + k 2 r + ki 



(C.2) 



for any natural base r. Those numbers are not important 
in themselves, what matters is thinking about the sets of 
{k m } as numbers in an abstract representation to base r. 

On entry to our algorithm we have an integer n > 
0, and on exit we wish to have the number nsol of all 
solutions of Eq. (|3l]) and the solutions themselves. For 
a set Si we introduce an abstract variable S to describe 
the array k ( 1 : n) of n integer elements ordered as "quasi- 
decimal digits" . 

So, in an abstract form the algorithm looks like: 



12 



S. Blinnikov & R. Moessner: Expansions for nearly Gaussian distributions 



Statement QS : 

all S <= S (CURRENT) 
are out if and only if 
k(l)+2*k(2)+. . .+n*k(n)=n 
<*initiate: make QS true for S=S (INITIAL) *>; 
_while <*condition: S *= S (FINAL) *> _do 
<*nextset: find next S 

keeping QS invariant *> 
_od; — Here S=S (FINAL) and QS=.TRUE., i.e. 
— all needed solutions are found 

It is easy to initiate QS: 

k(l)=n; 

mold=l; — keeps largest m 

— for non-zero k(m) 
_do m=2,n; 

k(m)=0 
_od; 
nsol=l ; 

The integer mold keeps track of the progress of the algo- 
rithm and allows us to establish the condition for the 
continuation of the main loop, 

mold 7^ n . 

Finally, the core of the algorithm is the node next set 
where we work with our sets as with digits, which is like 
doing the addition of decimal numbers by hand, column 
by column from right to left. 

advance S (CURRENT), i.e. consider adding 

n-( 2*k(2) +3*k(3) + . . . +mold*k(mold) ) 

to k(l) trying to add 1 to the lowest of 

k(2) , . . . ,k(mold) possible. 

If adding 1 to the lowest k(m) 

makes the sum > n then 

put this k(m)=0 and try k(m+l) 



Table C.l. Table of solutions to Eq. (|3 
Zeros for higher orders are left blank 



for n = 1 to 



n 1 


2 


3 


4 


1 


2 


3 


4 




10 


11 


12 






100 


20 








101 








1000 


n 5 


6 


7 


8 


5 


6 


7 


8 


13 


14 


15 


16 


21 


22 


23 


24 


102 


30 


31 


32 


110 


103 


104 


40 


1001 


111 


112 


105 


10000 


200 


120 


113 




1002 


201 


121 




1010 


1003 


202 




10001 


1011 


210 




100000 


1100 


1004 






10002 


1012 






10010 


1020 






100001 


1101 






1000000 


2000 








10003 








10011 








10100 








100002 








100010 








1000001 








10000000 



For reference, we present the first 8 sets of solutions 
found by this algorithm. In practice it is easier not to use 
the table below even for low n, but generate all coefficients 
in the code. 



integer sum of current set S 



m=l ; 

sumcur=n ; 
_repeat 

sumcur=sumcur-k(m) *m+m+l ; 

k(m)=0; 

k(m+l)=k(m+l)+l; 
m=m+l ; 

_until sumcur <= n _or m>mold; 
_if m>mold _then mold=m _fi; 
k (l)=n- sumcur ; 
nsol=nsol+l ; 

Here the node nextset ends, and the whole algorithm 
is finished. We have written it here in a pseudocode which 
can be translated mechanically into any machine language. 
We actually use a special preprocessor Tref or which au- 
tomatically transforms the text above to standard For- 
tran (see Weinstein & Blinnikov 1984, and Bartunov et 



al. |1997|) 
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